#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <sys/time.h>
#include <pthread.h>

/* 	Application of Monte Carlo simulation: Black Scholes model 	*/

long long  n_iteractions;		// number of iteractions for each thread
double * trials;			// vector of results

double S, E, r, sigma, T;		// parameters
long long M;

struct BoxMullerState
{
        double x1, x2, w, y1, y2;
        int useLast;
        struct drand48_data random;
};

void initBoxMullerState(struct BoxMullerState* state)
{
        state->random.__init = 0;
        state->useLast = 0;
        
        struct timeval now;
	gettimeofday(&now, NULL);
	state->random.__x[0] = now.tv_usec;
}

double boxMullerRandom(struct BoxMullerState* state)
{
        double randomNumber;

        if (state->useLast)
        {
                state->y1 = state->y2;
                state->useLast = 0;
        }
        else
        {
                do
                {
                        drand48_r(&state->random, &state->x1);
                        state->x1 = 2.0 * state->x1 - 1.0;
                        drand48_r(&state->random, &state->x2);
                        state->x2 = 2.0 * state->x2 - 1.0;
                        state->w = state->x1 * state->x1 + state->x2 * state->x2;
                }
                while (state->w >= 1.0);

                state->w = sqrt((-2.0 * log(state->w)) / state->w);
                state->y1 = state->x1 * state->w;
                state->y2 = state->x2 * state->w;
                state->useLast = 1;
        }

        return state->y1;
}

double max(double a, double b){
	if(a > b)
		return a;
	else
		return b;
}

double avg(){
	int i;
	double sum = 0;
	for(i = 0; i < M; i++)
		sum += trials[i];
	
	return sum/M;
}

double variance(double mean){
	int i;
	double sum = 0;
	for(i = 0; i < M; i++)
		sum += (trials[i] - mean) * (trials[i] - mean);
	
	return sum/M;
}

/*	Funcion used by all threads 	*/
void * function ( void * args )
{
	// initialize variables
	int * args_d = (int *) args;
	int thread_id = *args_d;
	long long i;
	double number, t;

	struct BoxMullerState state;

	initBoxMullerState(&state);

	for(i = n_iteractions*thread_id; i < n_iteractions*(thread_id+1); i++){
		
		number = boxMullerRandom(&state);

		t = S * exp( (r - (sigma * sigma) / 2.0 ) * T + sigma * sqrt(T) * number);

		trials[i] = exp( -r * T ) * max(t - E, 0);
	}

	return NULL;
}

int main (int argc, char ** argv)
{	
	if(argc != 2){
		printf("Usage: ./black_scholes-par <number_of_threads> \n");
		exit(1);
	}

	scanf("%lf %lf %lf %lf %lf %lld", &S, &E, &r, &sigma, &T, &M);

	/*	initialize variables	*/
	int n_threads = atoi(argv[1]);			// number of threads
	n_iteractions = M / n_threads;			// number of iteractions for each thread

	// vector of threads and thread ids 
	pthread_t * threads = (pthread_t *) malloc( sizeof( pthread_t) * (n_threads + 1)); 
	int * thread_id = (int *) malloc( sizeof( int ) * (n_threads + 1) );

	// vector of results
	trials = (double *) malloc(sizeof(double) * (M + 1) );

	// create all threads
	int j;
	for(j = 0; j < n_threads; j++) {
		thread_id [j] = j;
		if(pthread_create(&threads[j], NULL, function, &thread_id[j]))
			fprintf(stderr, "Error creating thread %d\n", j), exit(1);
	}

	// join all threads
	for(j = 0; j < n_threads; j++) {
		if(pthread_join(threads[j], NULL))
			fprintf(stderr, "Error joining thread %d\n", j), exit(1);
	}

	double mean = avg();
	double stddev = sqrt (variance(mean));
	double confwidth = 1.96 * stddev/sqrt(M);
	double confmin = mean - confwidth;
	double confmax = mean + confwidth;

	printf("S	%.lf\n", S);
	printf("E	%.lf\n", E);
	printf("r	%.lf\n", r);
	printf("sigma	%.lf\n", sigma);
	printf("T	%.lf\n", T);
	printf("M	%.lld\n", M);
	printf("Confidence interval: [%.6lf, %.6lf]\n", confmin, confmax);	

	// free memory
	free(threads);
	free(thread_id);
	free(trials);

	return 0;
}
